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'. Abstract 
(N . 

■ Molecular dynamics simulations of water, liquid beryllium fluoride and silica melt are 

^ I used to study the accuracy with which the entropy of ionic and molecular liquids can be 

d • estimated from atom-atom radial distribution function data. The pair correlation entropy is 

^ I demonstrated to be sufficiently accurate that the density-temperature regime of anomalous 

. behaviour as well as the strength of the entropy anomaly can be predicted reliably for both 

3 I ionic melts as well as different rigid-body pair potentials for water. Errors in the total 
thermodynamic entropy for ionic melts due to the pair correlation approximation are of the 

CN ■ 

'i> ' order of 10% or less for most state points but can be significantly larger in the anomalous 

^ ! 

G^ ■ regime at very low temperatures. In the case of water, the rigid-body constraints result 

■ 

I in larger errors in the pair correlation approximation, between 20 and 30%, for most state 

^ . points. Comparison of the excess entropy, Se, of ionic melts with the pair correlation entropy, 

^1 5*2, shows that the temperature dependence of Sg is well described by T^^/^ scaling across 

'i> ■ both the normal and anomalous regimes, unlike in the case of 5*2. The residual multiparticle 

^ I entropy, AS = Se — S2, shows a strong negative correlation with tetrahedral order in the 

. . anomalous regime. 



* Author for correspondence (Tel: (+) 91 11 2659 1510; Fax: (+) 91 11 2686 2122; E-mail: 
charusOchemistry . iitd . ernet . in) 

1 



I. INTRODUCTION 



The thermodynamic excess entropy of a fluid is defined as the difference in entropy 
between the fluid and the corresponding ideal gas under identical temperature and density 
conditions. In the case of a classical fluid, the excess entropy corresponds to the lowering 
of the entropy of the fluid, relative to the ideal gas, due to the presence of multi-particle 
positional correlationsi. The total entropy of a classical fluid can be written as^^dik^^dn^il^ 

oo 

S = S,d + S, + Ss + ... = S,d + J2Sn (1) 

ra=2 

where Sid is the entropy of the ideal gas reference state, Sn is the entropy contribution due 
to n-particle spatial correlations and the excess entropy is defined as, Se = S — Sid- Such a 
multiparticle expansion of the entropy is interesting because it allows for the prediction of 
thermodynamic properties from structural correlation functions. Moreover, semiquantitative 
excess-entropy based scaling relationships of the form 

X = Aexp{-aSe) (2) 

where X is a suitably scaled transport propertyii'i^'^^i'i^ii^ allow for the possibility of mak- 
ing interesting connections between structural, thermodynamic and transport properties of 
fluids. 

Neutron or X-ray scattering experiments can provide atom-atom radial distribution func- 
tions (RDFs), 5'a/3(r), associated with the probability of finding an atom of species /3 at a 
distance r from an atom of species a relative to the probability in the corresponding ideal 
gasi. In terms of the atom-atom radial distribution function, one can write an ensemble- 
invariant expression for the pair correlation contribution to the entropy of a multicomponent 
fluid of N particles enclosed in a volume V at temperature T as^^ 

POO 

S2/NkB = -27: p XgXp / {gap{r) In 5(^/3 (r) - [gap{r) - l\]r'^dr (3) 

where Xa is the mole fraction of component a in the mixture and p is the number density of 
the fluid. Note that the excess entropy is measured in this case with respect to the entropy 
of an ideal gas consisting of a non-interacting mixture of spherical particles with entropy Sid 
given by 

Sid/NkB = 1 + 5^ x„ [0.5t;„ - Hpa^D] (4) 



2 



where Va is the number of degrees of freedom associated with the component a and 
is the thermal de Broghe wavelength of component a. The multiparticle expansion of the 
entropy has been studied most extensively for simple liquids and liquid mixtures of struc- 
tureless particles with isotropic interactions^'^i^'^i^. For such systems, the pair correlation 
contribution captures 85-90% of the excess entropy, with the largest discrepancy occurring 
at intermediate densities^. 

In the case of a homogeneous, isotropic molecular fluids, an alternative and physically 
reasonable reference state for the multiparticle expansion is provided by an ideal gas of rigid 
rotors. In this case, pair correlation contribution to the entropy must be formulated in 
terms of orientationally-dependent pair correlation functions, g{r,uj'^), where cu^ denote the 
relative orientations of the two molecules^ii^. The orientation-dependent pair correlation It 
is convenient to perform a further decomposition of g{r,u'^) as 



where g{uj'^\r) is the conditional probability of observing an orientational separation uj"^ at 
a relative separation r of the molecular centres of mass. The pair correlation entropy of 
a molecular system, denoted by 5*™°', may then be written as a sum of the translational 
entropy, S2 , dependent only on g{r), and orientational contribution, involving a suitable 
weighted average over g{uj'^\r). While simulations show that this procedure yields reasonable 
estimates of the entropy for standard pair-additive water models, orientation-dependent 
pair correlation functions are not directly accessible from experiment. Since evaluation of 
the correlation functions from simulations with adequate statistical accuracy may prove 
complicated, a number of approximate factorizations of g{r,uj'^) have been studiecP^'^'^. 

Since atom-atom radial distribution functions (RDFs) are obtainable from experiments, 
obtaining estimates of thermodynamic and transport properties of liquids from atom-atom 
RDFs is very attractive because it allows one to directly connect structural information with 
thermodynamic and transport properties of fluids. In principle, such studies could contribute 
to improving reverse Monte Carlo strategies for structure determination as well as multi- 
scale methodologies for constructing coarse-grained potentials2^>2ii2^i2^i2Ii2^i2£. In order to 
develop this possibility, it is necessary to test the numerical validity of the pair correlation 
approximation for a wider range of systems than simple atomic liquids. The purpose of this 
paper is to examine the extent to which the pair correlation entropy, computed from atom- 
atom RDFs, captures the temperature and density dependence of the excess entropy for 
ionic melts and molecular liquids. We also study the behaviour of the residual multiparticle 




(5) 
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entropy (RMPE), defined as AS" = Se — S2, since it has been shown contain physically 
significant information on local order in a liquid close to phase transitions. For example, 
the state point for which the RMPE crosses from negative to positive is found to be close to 
ordering phase transitions, such as freezing and the isotropic- nematic transitiort^*^'^'^. 

Ionic melts are of considerable technological importance and are, as a consequence, well- 
studied experimentally as well as computationally^^i^. Many ionic melts form random liquid 
state networks and show interesting deviations from simple liquid behaviour in their ther- 
modynamic and transport properties2ii2^i2^i21i2^. Other than a limited comparison in the 
case of liquid silica^^i^i^, we are not aware of any systematic testing of the pair correlation 
approximation for the entropy in the case of ionic melts. As representative examples of ionic 
melts, we choose beryllium fiuoride (BeF2) and silica (Si02). Molecular liquids represent 
an interesting test case for estimating entropy using atom-atom RDF data, since the pair 
correlation approximation, as defined in equation (3), is expected to be poorer because the 
bond angle contraints give rise to strong three-body correlations and an ideal gas of rigid 
molecules is a more appropriate reference state, than a mixture of monoatomic ideal gases 
as implicit in equation (3). Nonetheless, the experimentally accessible quantities are the 
atom-atom RDFs and it is interesting to see how much information can be extracted from 
these quantities. As an example of a molecular liquid, we choose water which is of obvi- 
ous interest as a pure liquid as well as a solvent, and is representative of a wider class of 
hydrogen-bonded solvents^i^i^i^. Computing the entropy of a molecular fiuid from atom- 
atom pair correlation data has not been attempted, though orient at ionally dependent pair 
correlation functions, which are not accessible experimentally, have been used to estimate 
the entropy of several commonly used model water potentials^id^. 

All three liquids that we study in this work (H2O, Si02 and BeF2) have similar anomalous 
thermodynamic and kinetic properties^ii^iS^i^iSSii^i^didld^dSi^iMii^i^. The most obvious 
signature of the thermodynamic anomalies is the existence of a regime of anomalous density 
behaviour where the isothermal expansion coefficient, a = {l/Vm){dVm/dT)p is negative. 
The region of the density anomaly is bounded by temperatures of maximum and minimum 
density for which a = 0. While the temperature of minimum density is experimentally 
difficult to observe, the locus of temperatures of maximum density (TMD) can be traced in 
the density-temperature or temperature-pressure plane for many systems. The maximum 
temperature and maximum density along this locus, and p^md respectively, can be 

thought of as upper bounds above which the thermal kinetic energy or degree of compression 
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respectively is sufficient that the hquid conforms to thermodynamic behaviour characteris- 
tic of simple liquids. Note that the lowest temperature or maximum density that can be 
obtained on the TMD locus will be limited by spontaneous crystallisation or glass forma- 
tion. The diffusional anomaly corresponds to the set of state points for which the diffusivity 
increases with density. A unified explanation for the thermodynamic and kinetic anomalies 
can be found based by existence of an excess entropy anomaly involving a rise in excess 
entropy on isothermal compression^^i^i^i^i^i^. The relationship between excess and pair 
correlation entropies in these liquids is of special interest from the point of view of comparing 
the two quantities in the normal and anomalous regimes, since there is some indication that 
anomalous colloidal fluids with core-softened interactions show a greater magnitude of the 
RMPE contribution at high densities^. In the case of water, simulation results for TIP4P 
water suggest that the RMPE is close to zero close to the TMD at 1 atm pressure. 

Molecular dynamics simulations for BeF2 and Si02 are performed using the transferable 
rigid ion model (TRIM)^i^i^i^ and van Beest-Kramer-van Santen (BKS}^i^i^ potentials 
respectively. For both the systems, the ideal gas limit is defined as an MX2 binary mixture 
of uncharged particles with the same masses as the corresponding ionic melts but with zero 
interparticle interactions. The excess entropy, Se, with respect to this ideal gas is evaluated 
using thermodynamic integration. Using the radial distribution functions extracted from 
the simulations, one can obtain both the pair correlation entropy (5*2) and the residual 
multiparticle entropy (RMPE), given by AS = — 82- In the case of water, we have 
simulated three of the commonly used rigid-body, effective pair interaction models(SPC/13^, 
TIP3P^^ and TIP5P^^) and evaluated the entropy from the gooi.^)^ 9HH{.f) and gonir) radial 
distribution functions. Table I shows the and Ptmd values for BKS silica, TRIM BeF2 

and the SPC/E and TIP5P water models as determined in previous studies^i^i^i^i^. 

The paper is organised as folllows. Section II summarises the computational details 
associated with the molecular dynamics simulations of silica, beryllium fluoride and water 
that we have performed. Section III describes the thermodynamic integration procedure 
that we have followed to obtain the thermodynamic entropy of the two MX2 ionic melts. 
In the case of water, thermodynamic entropies for standard potentials are available in the 
literature. Section IV discusses the results for the two ionic melts while Section V contains 
the results for water. Conclusions are summarised in Section VI. 
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II. COMPUTATIONAL DETAILS 



For all the three systems, we use parametric effective pair potentials. Each of these inter- 
action models has been well-studied in the literature in terms of behaviour on supercooling, 
glass transition and waterlike structural, kinetic and thermodynamic anomalies. For con- 
venience, we describe the functional forms of the pair potentials and give the parameters 
in Tables II, III and IV. The molecular dynamics simulations were performed using the 
DL_POLY software packag 



1. Model Potential for Silica 

To model the interatomic interactions in silica, we use the van Beest-Kramer-van Santen 
(BKS) potential with an additional 30-6 Lennard- Jones type correction term^i^i^. The 
pair interaction between atoms i and j is given by: 

30 



>BKs[r,,) ^^^^^^^ + A,,exp + 4e,,, 



(Tij 



(6) 



where is the distance between atoms i and j carrying charges and g^, A^-, hij and 
Cij are the parameters associated with the Buckingham potential for short-range repulsion- 
dispersion interactions and and oij are the energy and length scale parameters for the 
30-6 Lennard- Jones interaction. The parameters for the modified BKS potential used in this 
work are given in Table II. 



2. Model Potential for Beryllium fluoride 

We use the transferable rigid ion model (TRIM) potential for interatomic interactions in 
beryllium fluoride^*^'^'^. The pair interaction between atoms i and j is given by: 

0™(r.,) = + f 1 + ^ + ^"j 6exp (^.l±ll^Il2\ (7) 

where r^j is the distance between atoms i and j. The parameters associated with an atom of 
type / are the charge zi , the number of valence-shell electrons rii and the ionic size parameter 
ai- The repulsion parameter h and the softness parameter s are assumed to be the same for 
all three types of pair interactions in BeF2. The parameters for the TRIM potential used in 
this work are given in Table III. 
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3. Model Potentials for Water 



All rigid-body, effective pair potentials for water assume that the molecule can be rep- 
resented by a single Lennard- Jones site located on the oxygen atoms and model the charge 
distribution by a set of distributed charges with a fixed geometry. The parametric form of 
the interaction between two water molecules a and b is given by: 



where i and j index partial charges located on molecules a and b respectively and roo refers 
to the distance between the oxygen atoms of the two monomers. 

The SPC/E and TIP3P models use three charged sites which are located at the atomic 
positions and are associated with the corresponding atomic masses. The TIP5P model 
retains the three atomic sites and uses two additional massless sites to represent the oxygen 
lone pair electron density distribution. The potential parameters for the three water models 
are summarised in Table IV. 

4. Molecular Dynamics 

Molecular Dynamics simulations for all three systems were carried out in the canonical 
(N-V-T) ensemble, using the DL_POLY software package^>^, under cubic periodic boundary 
conditions. The effects of electrostatic (long-range) interactions were accounted for by the 
Ewald summation method^"^. The non-coulombic part was truncated and shifted at 7. 5 A 
for Si02 and BeF2 and at 9.0^4 for H2O. A Berendsen thermostat, was used to maintain 
the desired temperature for the production run. The leapfrog Verlet algorithm with a 
time step of Ifs was used to integrate the equations of motion; in the case of water, the 
SHAKE algorithm was used to impose the rigid-body contraints. The MD simulation details, 
including the lengths of the equilibration and production runs, are summarised in Table V. 
The cutoff distance for the radial distribution functions was kept the same as the potential 
cutoff and the bin size was adjusted so that trapezoidal qaudrature to evaluate the pair 
correlation contribution to the entropy using equation (3) resulted in less than 1% error. 
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III. ESTIMATION OF THERMODYNAMIC ENTROPY OF IONIC MELTS 

We follow the thermodynamic integration procedure described previously by Saikia- 
Voivod and co-workers to evaluate the thermodynamic excess entropy of BKS silica 
where the full Ewald summation was replaced by smoothly tapering off to zero the real- 
space and short-range contributions to the interaction and removing the reciprocal space 
contribution^i^i^i^. The internal energy differences between the short-range and full Ewald 
versions of the BKS were of the order of 0.2% or less. Since we use the full Ewald version 
of the potential, it was, however, necessary to recompute the thermodynamic entropy of 
BKS silica. To our knowledge, the entropy of the TRIM model of BeF2 has so far not been 
evaluated. 

The model potentials used here envisage the ionic melts as being fluids composed of 
particles with long-range Coulomb interactions and short-range repulsions. Such a Coulom- 
bic system cannot be reversibly transformed into an ideal gas reference system. Therefore, 
a binary Lennnard- Jones (BLJ) system with an MX2 stoichiometry was introduced as an 
intermediate state. The entropy of the ionic melt at a given state point (Vo,To) was first 
evaluated relative to the binary Lennard- Jones system at the same state point. The entropy 
of the BLJ liquid at (Vo,To) was then evaluated relative to the ideal gas limit at 
where V^o corresponds to such a large volume that the fluid can be regarded as an ideal gas. 

To determine the entropy of the ionic melt relative to the BLJ liquid, thermodynamic 
integration was performed with respect to a Kirkwood-type coupling parameter. A, which 
interpolates between the potential energy function, f/(r), of the BLJ liquid and the ionic 
melt such that^- 



The free energy difference between the BLJ and MX2 ionic melt at {Vq,Tq) can be written 
as: 



The integrand (■ ■ ■)a in the above equation refers to a canonical ensemble average for a 
system described by the potential energy function Ux{r). The integral was evaluated using 
trapezoidal quadrature with the value of the integrand at specific values of the variable A 
being computed using a canonical ensemble MD simulation with volume Vq, temperature Tq 
and interaction potential Ux{r). The entropy of system MX2 at this state point, denoted 




(9) 




(10) 
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by Smx2{Vo,To), then be evaluated using, 

^MX,(^o,ro) = Sblj{Vo,To) + i [Umx,{Vo,To) - Ublj{Vo,To) - AF{Vo,To)] (11) 
The entropy of the BLJ system relative to the binary ideal gas is obtained as 



1 rVoc 
SBLj{Vo,n) = Sid{Vo,To) + - 



{Ublj{Vo,To)) - PIW 



(12) 



'Vo 

where Sid is the entropy of the ideal gas as defined in equation (4) The third term in the above 
expression is obtained by evaluating the excess pressure of the BLJ system for increasing 
system volumes using simulations at moderate to low densities and a virial expansion at 
very low densities. 

Once the entropy of the ionic melt has been obtained at a given state point (Vq, Tq) using 
the above procedure, the entropy at an arbitrary state point (V, T) can be obtained as 

S{V, T) = SiVo, To) + A5r + A5y (13) 

where /S.St is entropy change for isothermal change in volume from Vq\joV and /S.Sv is the 
entropy change for isochoric change in temperature from Tq to T. The following expression 
can be used to evaluate A5't: 

-V 



To 



Umx,{V,To)-UmxM,To)+ [ P{V')dV' 



'Vo 

where P is the pressure of the MX2 hquid. To evaluate A5'y, we use 



(14) 



where E is the total internal energy which must be the sum of the classical kinetic (l.bNkBT) 
and potential {{UMX2)) energy contributions for the N particle system. Since the ensemble 
average of the potential energy for BKS silica as well as TRIM BeF2 can be fitted by 
functional form, a{V) + 6(V)T^/^, the integrand in equation (14) is simple to evaluate. 

Since the internal energy estimates obtained using the short-range BKS and the BKS 
potential with explicit Ewald summation used here are very small, we used the same BLJ 
reference state as used by Saika-Voivod and co-workers. The potential energy parameters 
for this BLJ system, referred to as BLJ5j02) listed in Table VI. and the entropy of this 
system at 4000K and 2.307 g cm~^ was previously evaluated as 84.028 J mol~^ K~^. For 
this state point, we used equations (13) and (14) in conjunction with molecular dynamics 
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simulations of systems at different values of A to compute the entropy of BKS silica as 76.82 
J K~^, which differs by just 2.5% from the entropy of short-range BKS sihca. The entropies 
of BKS silica over the entire range of temperature and density were then calculated relative 
to this state point. 

Given the stoichiometry and the similar radius ratios of BeF2 and Si02, it was possible to 
rcscale the BLJsi02 order to define an appropriate reference system for BeF2, denoted by 
BLJBeF2- Since the locations of the first minima in the Be-Be, F-F and Be-F RDFs are very 
similar to those in the corresponding RDFs for silica, the length scale parameters for the 
BL J BeF2 state system were kept the same as for BL J5i02 ■ The relative well-depths of the three 
LJ-type pair interactions were also kept the same. Therefore if cm-x and aM-x are chosen 
as the units of energy and length respectively, then the two BLJ systems would have identical 
internal energies and entropies at the same number density and reduced temperature. To set 
a suitable energy scale for the BLJ system, we use ese-F ~ e5i-o(r™l),seF2)/(^TMX),Si02) 
and set e^e-F as 14.21K, as shown in Table VI. As discussed above, the BLJ5j02 liquid 
at 4000K and 2.307 g cm"''^ has total entropy equal to 84.028 kJ mol"^ K~^= lO.lOeA;^, 
excess entropy of 9.23A:^ and internal energy of 56.40e5j_o- The equivalent state point for 
BeF2 melt would be at 1777K and 1.805 g cm^^ which would have the same internal energy 
and excess entropy in reduced units. In order to calculate the total entropy, the ideal gas 
contribution must be computed using masses of Be and F. The total entropy of BLJ bcF^ was 
then found to be 70.7 J mol~^K~^=8.5/cB at this state point. The entropies over the entire 
range of temperature and density for BeF2 were then calculated relative to this state point. 



IV. IONIC MELTS: SILICA AND BERYLLIUM FLUORIDE 

Figure 1 shows the thermodynamic excess entropy as a function of density for Si02 and 
BcF2 melts. Unlike the monotonic decrease in with p seen in simple liquids, both systems 
have a well-defined anomalous regime where excess entropy rises with increasing density. 
The anomalous behaviour is more pronounced for the low temperature isotherms where the 
anisotropic interactions stabilizing local tetrahedral order are large in comparison to thermal 
kinetic energies. All the Se{p) curves have a well-defined maximum beyond which the excess 
entropy decreases with density; the maxima are more pronounced for the low temperature 
isotherms with T < T^md but can be identified in the high temperature isotherms as well. 

It is useful to compare the Se{p) behaviour with the density dependence of the pair 
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correlation entropy, 5*2, for Si0222 and BeF2^. The maxima in Se{p) and 5*2 (p) curves 
are located at essentially the same value of density within simulation error. For example, 
maxima in Se coincides with the maxima observed in 5*2 at p = 3.4g cm~'^ at 4000 K for Si02 
and at p = 2.82 g cm~^ at 1500 K for BeF2. The most notable qualitative difference between 
the Se{p) and 5*2 (p) curves is the absence of a minimum in the Se curves at low densities. 
The 5*2 (p) curves for all waterlike systems show a minimum at low densities which coincides 
with the maximum in the tetrahedral order parameter along the isotherms, indicating that 
5*2 is more sensitive to tetrahedral order than Se- 

Figure 2 shows the temperature dependence of the excess entropy of BeF2 and Si02 
along different isochores. We plot Se as a function of T~^/^ since earlier density functional 
as well as simulation studies suggest that the excess entropy of simple liquids, even within 
the pair correlation approximation, obeys a T'"^/^ scaling^*^. Figures 2(a) and (b) show 
that both Si02 and BeF2 melts show an approximate T^^'^^^^ scaling, with small deviations 
from the scaling behaviour in the anomalous regime. Note that a T^^^ scaling predicted 
for the configurational potential energy is obeyed by both the melts and has been used 
when performing thermodynamic integration^^. Interestingly, the pair correlation entropy, 
5*2, plotted as a function of T~^/^, of liquids with waterlike anomalies shows a significant 
change between on going from the normal to the structurally anomalous regime. 

To summarise, the qualitative differences that emerge in the density and temperature de- 
pendence of Se and 5*2 are: (i) absence of a minimum in Se, unlike the minimum in 5*2 which 
correlates with the maximum in tetrahedral order and (ii) the temperature dependence of 
Se along isochores is well described by T^^/^ scaling across both the normal and anomalous 
regimes, unlike in the case of 5*2. This suggests that it is important to examine the relation- 
ship between the residual multiparticle entropy, AS, and tetrahedral order parameter, qtet, 
as is done below. 

Figure 3(a) displays the correlation between pair correlation entropy 5*2 and thermody- 
namic excess entropy, Se, in the case of of Si02. A strong correlation between 5*2 and Se 
is seen for the low isotherms of T = 5000, 4500 and 4000 K i.e. for T/TJ7^ < 1. For 
temperatures greater than (T = 5500 and 6000 K), the 5*2 versus Se curve shows two 

branches, a high density and a low density one. The density demarcating the two branches 
corresponds to 3.0 g cm^^. Figure 3(b) shows the 5*2 versus Se correlation plot for BeF2. 
Isotherms above 2000K {T^^fj-, = 2310/^), the curves have distinct high and low density 
branches with the demarcation density being 2.4 g cm~^. 
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In order to understand the relationship between Se and 5*2, we plot AS* as a function of p 
in Figure 4. Several interesting features of the RMPE as a function of density are common 
to both BeF2 and Si02 melts. In both systems, AS" is essentially constant after a charac- 
teristic density Ptmd which marks the maximum density for which the thermodynamically 
anomalous behaviour in the density is observed in the liquid phase. This corresponds to a 
density of 3 g cm~'^ for Si02 and 2.4 g cm~^ for BeF2. In both systems isotherms which lie 
at or below show a decrease in AS* with increasing p till Ptmd is reached. In contrast, 

isotherms which lie above T^^f^ show an increase in AS till p^md reached. 

Figure 4 also shows that the RMPE is negative for all state points in the case of Si02 
but not in the case of BeF2. This may in part be due to the fact that the BKS interaction 
potential in the case of Si02 does not include any explicit cation-cation (or Si-Si) interactions. 
We find that the RMPE contributes about 10% or less to the excess entropy, Se, of the ionic 
liquids in the normal regime. In the anomalous regime, specially at low temperatures, the 
contribution of the AS terms can be as large as 30% of the total excess entropy in terms of 
the absolute magnitudes. 

To assess the quantitative reliability of the pair correlation approximation to obtain the 
total thermodynamic entropy. Figure 5 shows the quantity |AS'|/S'% where S = Sid + Se is 
the total thermodynamic entropy. The error caused by the pair correlation approximation 
is of the order of 6% or less for both ionic melts in the normal regime, but as expected from 
the behaviour of the RMPE, the error can be greater in the anomalous regime. Other than 
the 1500K isotherm of BeF2, the magnitude of the error for all state points studied is less 
than 10%. 

The above results suggest that the nature of the relationship between the residual multi- 
particle entropy and structural order is significantly different in the normal and anomalous 
regime. Therefore, we plot AS" as a function of the tetrahedral order parameter qtet in Figure 
6. The local tetrahedral order parameter, qtet, associated with an atom i is defined as 



where ifjjk is the angle between the bond vectors rij and where j and k label the four 
nearest neighbour atoms of the same type^. Figure 6 shows that AS has a strong negative 
correlation with tetrahedral order in this regime. Clearly, the dominant contribution to 
the RMPE in the anomalous regime is from the three-body terms which are strongly anti- 
correlated with tetrahedral order. 




(16) 



j=i k=j+i 
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We now reconsider the behaviour of 5*2 versus Se at temperatures less than or equal to 
"^TMD Figure 3. 5*2 and Se are strongly correlated and increase with increasing density till 
Pmax is reached. The anomalous regime also shows a strong correlation between tetrahedral 
order and 5*2 or any related measure of structure in the pair correlation function^. The 
effect of the local tetrahedral ordering, imposed in the case of ionic melts from the steric 
factors associated with the relative ionic radii, is to strongly couple the two and three-body 
contributions to the excess entropy. For isotherms lying above T^"^, at low densities 5*2 
is essentially constant while Se increases. The multiparticle correlations in this case serve 
to increase the entropy with increasing density and are essentially uncorrelated with the 
tetrahedral order. For densities above pmax and temperatures above Ti^^j-,, the system 
behaves essentially as a simple hquid, with 5*2 and Se both showing a negative correlation 
with density though with a relatively small variation in either entropy or local order. 



V. MODEL POTENTIALS FOR WATER 



We first consider the the pair correlation estimate of the entropy, based on the goo{f^), 
gonif') and qhh pair correlation functions, with the total entropy as estimated from thermo- 
dynamic integration for the SPC/E model for water. The experimentally obtained goH{f) 
and gHH{r) RDFs will show a sharp, narrow peak corresponding to the intramolecular dis- 
tances between these atoms. In the case of rigid body molecular dynamics, this peak will 
reduce to a 5-function. The 5*2 contribution can be computed from equation (3) from the 
g{r) functions without this additional 5-function contribution. 

Figure 7(a) shows the density dependence of the thermodynamic entropy, S*, taken from 
ref.— , for several isotherms of SPC/E water lying between 210K and 300K. Comparison 
with the pair correlation estimate, S* ~ Sid + 'S'2, shown in Figure 7(b), indicates that the 
pair correlation entropy shows the correct qualitative behaviour. The anomalous rise in 
entropy between 0.95 and 1.10 g cm~^ is seen in both Figures 7(a) and 7(b) for isotherms at 
210,220 and 230 K. The quantitative errors introduced by the pair correlation approximation 
for SPC/E water are assessed by plotting \ AS\/S% where S is the thermodynamic excess 
entropy taken from ref.— as a function of density in Figure 7(c). The errors due to the 
pair correlation approximation at temperatures of 240K and above lie between 20 and 30% 
and are almost independent of density. Not surprisingly, this is larger than errors observed 
for the ionic melts in the normal liquid regime. At low temperatures, specially in the 
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anomalous regime, there is a strong density dependence of the errors; for example, for the 
21 OK isotherm, the error varies between 10 and 60% . 

Since the pair correlation entropy provides a very reasonable estimate of the temperature- 
density regime in which thermodynamic anomalies may be expected, we consider the two 
other effective pair potential models of water known to have very different regimes of anoma- 
lous behaviour (Table I). The TIP5P model reproduces the experimental data closely in 
this respect and shows a value of 282K. In the case of the SPC/E model, however, 

the anomalous regime is shifted to temperatures approximately 30 to 40K below those see 
experimentally. In the case of TIP3P, the anomalies appear to occur at even lower temper- 
atures though a detailed mapping of waterlike anomalies has not been performed for this 
model. The TIP3P model is, however, of considerable interest since it is frequently used 
to model the aqueous solvent in biomolecular simulationsSS. To illustrate the quantitative 
differences between the models, one may note that a recent study indicates that the TMD 
at 1 atm pressure occurs at 241K, 182K and 285K for the SPC/E, TIP3P and TIP5P mod- 
els respectively^. Figure 8 shows the density dependence of the pair corrrelation entropy, 
5*2, as a function of density for different isotherms for all three models. It is immediately 
apparent that the thermodynamically anomalous regime can be correctly identified for all 
three models from this pair correlation information. An excess entropy anomaly exists in the 
TIP3P, SPC/E and TIP5P models below 200K, 260K and 360K respectively. The strength 
of the excess entropy anomaly, as identified by {dS2/d{ln p))j^^ , is maximal in the TIP5P 
model and least in the TIP3P model. 



VI. CONCLUSIONS 



This paper compares the pair correlation entropy (52), determined from the atom-atom 
radial distribution functions with the excess entropy of two ionic melts (liquid silica and 
beryllium fluoride) and a molecular liquid (water). The three liquids that we have chosen to 
study show distinct normal and anomalous regimes. The anomalous regime is characterised 
by significant departures from simple liquid behaviour and shows a rise in thermodynamic 
entropy on isothermal compression. This anomalous entropy behaviour can be connected 
to the existence of waterlike structural, kinetic and thermodynamic entropies. The pair 
correlation entropy is sufficiently accurate that the density-temperature regime of anomalous 
behaviour as well as the strength of the entropy anomaly can be predicted reliably for both 
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ionic melts as well as different rigid-body pair potentials for water. This simple connection 
between atom-atom radial distribution functions and thermodynamic anomalies predicted 
by different water models has not been discussed previously in the literature. 

To assess the quantitative contribution of the pair correlations to the thermodynamic 
entropy, we compare the RDF-based 5*2 estimator for ionic melts with the thermodynamic 
excess entropy, Se, measured with respect to an ideal gas mixture. Errors in the total 
thermodynamic entropy for ionic melts due to the pair correlation approximation are of the 
order of 6% in the normal regime but can be significantly larger in the anomalous regime. In 
the case of water, we compared the total thermodynamic entropy with the estimate based on 
the atom-atom RDFs. As expected given the rigid-body constraints for molecular liquids, 
the pair correlation approximation then causes significantly larger errors, between 20 and 
30%, in the normal liquid regime. 

In the case of ionic melts, the high density limit of the anomalous regime is marked by 
the maximum in the Se{p) curves at a given temperature and this is well reproduced by the 
5*2 (p) curves. The excess entropy curve, Se{p), does not show a minimum, unlike the 5*2 (p) 
curves for which the minimum correlates with the maximum in tetrahedral order. Along an 
isochore, the temperature dependence of Se is well described by T~^/^ scaling across both 
the normal and anomalous regimes. In contrast, 5*2 shows T~^/^ scaling behaviour only in 
the normal liquid regime. 

Comparison of pair correlation approximation to the excess entropy from scattering data 
with calorimetric estimates of the entropy can provide interesting insights into the role of 
multiparticle interactions in different liquids. For example, in the case of ionic melts, we 
show that the relationship between 5*2 and Se is significantly different in the normal and 
anomalous regime. Strong, local anisotropic interactions are dominant in the anomalous 
regime and serve to couple the two- and three-body contributions to the entropy. Thus 
the correlation between local order metrics and RMPE in simple and anomalous liquids is 
qualitatively different. 

It is useful at this point to compare our results for excess entropy of the three tetrahedral 
network- forming liquids with the earlier results for H2O using an ideal gas of rigid rotors 
as a reference stai^^^. Clearly the separation of time scales between intramolecular and 
intermolecular interactions of water justifies the rigid-body ideal gas limit as the physically 
more reasonable choice. In the case of ionic melts, the molecular limit is not an obvious 
choice and defining the reference state as a non-interacting mixture of atoms is a reasonable 



15 



option. The greater accuracy of the pair correlation estimator for the excess entropy of ionic 
melts, as compared to water, is a consequence of the appropriateness of choice of reference 
state. Both approaches indicate that the contribution of multiparticle interactions to the 
excess entropy changes qualitatively on going from the normal to the anomalous regime. The 
translational entropy, as defined in refs.— corresponds to the contribution of the 0-0 pair 
correlation function to the pair correlation entropy as defined in our work. The orientational 
entropy must correspond to a subset of the three-body terms in our formulation which is 
supported by the strong correlation of the RMPE with the tetrahedral order parameter. 

In conclusion, the results presented in this paper suggest that atom-atom radial distri- 
bution functions can be used to construct a reliable, semiquantitative structural estimator 
for the entropy. The comparison of such a structural estimator for the entropy with calori- 
metric data could lead to interesting insights into the role of pair and higher-order particle 
correlations in determining thermodynamic and transport properties of liquids and serve 
as additional inputs for improving structure prediction by reverse Monte Carlo or related 
techniques based on RDF data. We also note that there has been considerable interest in 
developing coarse-graining strategies based on effective pair potentials for complex liquids. 
For example, isotropic models of water which reproduce the gooi'f') function have been con- 
structed but do have some problems with representability and transferability^i2^i21i2^i2£. It 
may be interesting to consider coarse graining strategies which modify the RDF data to 
reproduce the excess entropy of the liquid. 
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TABLE I: Maximum temperature and corresponding density for locus of state points corresponding 
to temperatures of maximum density (TMD) for different isobars. These temperatures and densities 
mark the maximum in the density-temperature (p — T)-plane of the curve connecting state points 
with {dV/dT)p = 0. The data for and P^md were taken from the hterature and the 

apprpriate reference is cited in the column heading. 

SPC/E53 TIPSP^'^ BeF2^ SiO^ 
T:p^f^/K 249 282 2310 4940 

Ptmd/ g cm-3 0.97 1.0 1.805 2.307 



TABLE II: Potential parameters for BKS Si02 with 30-6 Lennard-Jones correction term o i 



i j Aij bij Cij €ij (Tij 

(kJ mol-i) (i-i) (kJ mol-i) (kJ mor^) (i) 



0-0 134015 2.76 16887.3 0.101425 1.7792 
Si-0 1737340 4.87 12886.3 0.298949 1.3136 



TABLE III: TRIM potential parameters for 
p b 
(i) (i) (i-i) (kJmol-i) 



0.93 1.24 0.29 2 -1 2 8 34.33 
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TABLE IV: Comparison of parameters of common water modela^i^i^. SPC/E and TIP3P are 
three site models with the respective charges and masses centered at the O and H sitesS..13^ TIP5P 
is a five site model with two additional massless sites M^. 





SPC/E 


TIP3P 


TIP5P 


e (kcal mol ^) 


0.155 


0.152 


0.160 


a{A) 


3.166 


3.150 


3.120 


roH 


1.000 


0.9572 


0.9572 


ZHOH{deg) 


109.47 


104.52 


104.52 


QO 


-0.8472 


-0.834 


0.0 


QH 


0.4238 


0.417 


0.241 


QM 






-0.241 


roM{A) 






0.70 


ZMOM{deg) 






109.47 
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TABLE V: Computational details of molecular dynamics simulations of Si02, BeF2 and H2O. The 
simulation cell size is given in terms of the number of formula units, M, present in the system. 
An MD time step of 1 fs was used for all three systems. The time constant for the Berendsen 
thermostat is denoted by tb- Equilibriation and production run lengths are denoted by teq and 
tprod respectively. 



System 


Model M 




tprod 




MD algorithm 






(ns) 


(ns) 


(ps) 




Si02 


BKS 150 


3-6 


5-10 


200 


Verlet 


BeF2 


TRIM 150 


4-8 


6-8 


200 


Verlet 


H2O 


SPC/E 256 


1 


1 


10 


Verlet + Quaternion 


H2O 


TIP3P 256 0.25-0.5 0.5-0.75 


1 


Verlet + SHAKE 


H2O 


TIP5P 256 0.25-0.5 0.5-0.75 


1 


Verlet + Quaternion 



TABLE VI: Potential energy parameters for the binary Lennard-Jones fluid used as reference state 
for thermodynamic integration to obtain entropies of BKS Si02 and TRIM BeF2. 





SiOs 






BeF2 




i- 3 


e 


ad 


i- 3 


e 


a 




(kJ mol-i) 


(i) 




(kJ mol-i) 


{A) 


Si-0 


32.0 


1.6 


Be-F 


14.21 


1.6 


Si-Si 


23.0 


3.3 


Be-Be 


10.22 


3.3 


0-0 


23.0 


2.8 


F-F 


10.22 


2.8 



Figure Captions 

1. Thermodynamic excess entropy, 5*6, as a function of density for (a) Si02 and (b) 
BeF2. The vertical line shows the position of minima in the 5*2 (p) curves for SiO^^ 
and BeF2^. Unless otherwise stated, entropy is reported in units of A;^ per atom in 
all the figures. Isotherms are labelled in degrees K in all the figures. 

2. Temperature dependence of the thermodynamic excess entropy, S'e, for (a) Si02 and 
(b) BeF2 along different isochores. The isochores are labelled by the density in g cm~^. 
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3. Correlation between the pair correlation entopy, 5*2, and the thermodynamic excess 
entropy, Se, for (a) Si02 and (b) BeF2. The vertical arrows indicate the lowest density 
for each isotherm studied in this work. 

4. Variations in residual multiparticle entropy, AS, with density, p, along different 
isotherms, labelled by temperature in Kelvin, for (a) Si02 and (b) BeF2. AS is 
reported in units of ks per atom. 

5. Percentage contribution of residual multiparticle entropy to the total thermodynamic 
entropy, \AS\/S%, along different isotherms for (a) Si02 and (b) BeF2. 

6. Correlation between the tetrahedral order parameter, q and the residual multiparticle 
entropy, AS*, for (a) Si02 and (b) BeF2. AS is reported in units of ks atom. The 
arrows indicate the lowest density state point along an isotherm simulated in this 
study. 

7. Total entropy of SPC/E water as a function of density, (a) Total thermodynamic 
entropy, 5* taken from ref.^^ and (b) total pair correlation entropy, S* = Sjd + >S'2 
calculated by considering water as a binary mixture of H and O atoms. From top to 
bottom, the isotherms correspond to T = 300(v), 260(A), 240(A), 230(«), 220(0) 
and 210 K(B). Entropies are reported in units of ks per atom. 

8. Pair correlation entropy as a function of density for different model potentials of water: 
(a) SPC/E (b) TIP3P and (c) TIP5P. Entropies are reported in units of ks per atom. 
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FIG. 2: 
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FIG. 3: 
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FIG. 4: 
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FIG. 5: 
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FIG. 6: 
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FIG. 8: 



Figure 8 (a) 




(a) SPC/E 
21 OK T 
220K --7- 
240K * 
y 260K A 
■ 280K 

300K -o- 
* 320K ■ 
^ 340K B 



0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 



P (g cm"'') 



Figure 8 (b) 



l:ri 




P (g cm"'') 



Figure 8 (c) 
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